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Abstract 

We develop an approximate theoretical method to study discrete stochastic birth and death 
models that include a delay time. We analyze the effect of the delay in the fluctuations of the 
system and obtain that it can qualitatively alter them. We also study the effect of distributed 
delay. We apply the method to a protein-dynamics model that explicitly includes transcription 
and translation delays. The theoretical model allows us to understand in a general way the interplay 
between stochasticity and delay. 

PACS numbers: 
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I. INTRODUCTION 

Fluctuations play an important role in many areas of science [lj], and their study has 
become a well defined discipline. Delay in the interactions is also a common phenomenon 
in natural and artificial systems, and it is well known that it can alter qualitatively the 
dynamical behavior, for example inducing oscillations or even chaos [2j. In particular, both 
fluctuations and delay are relevant in gene-regulation systems, where a lot of research effort, 
both theoretical and experimental, has been recently carried out [2 2 1 . 

The connived effect of stochasticity and delay is not completely understood. In this 



context, many studies have 



Fokker- Planck equations 



time 



0,3 



ocused in delayed stochastic, Langevin, differential equations or 



that assume continuous variables, or random walks in discrete 



10] . Stochastic models with continuous time but discrete variables are the natural 



description of many systems such as chemical reactions, population dynamics, epidemics, 



etc. In some cases this discreteness is a major source of fluctuations 

In this work, we study some general stochastic birth and death processes that include 
delay. We follow a master equation approach that considers discrete variables in continuous 
time. We present an analytical treatment that allows us to study the effect of delay and show 
that the delay can alter qualitatively the character of the fluctuations. We also consider 
the situation with distributed delay and study how the fluctuations change as the delay 
distribution becomes wider. 

The paper is organized as follows: In the following section [III we present the theoretical 
approach, applying it to a general one-step birth-death model and discuss the influence of 
the delay in the fluctuations of this system. In section HTTl we consider the effect of distributed 
delay. In section IIVI we study a two-step transcription-translation model, more relevant to 
gene regulation. We finish in section [V] with some conclusions and comments. 

II. STOCHASTIC CREATION WITH DELAY 

Let us start by considering a simple one-step stochastic process in which the number of 
units (e.g. molecules) n of some compound X can only increase or decrease by one: 

C D 

0^X, X^0. (1) 
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The creation, C(n), and annihilation, D(n), rates depend, in general, on n. The probability 
P(n,t) that there are n molecules at time t follows a master equation[l|: 

= (E — l)[D(n)P(n,t)} + (ET 1 - l)[C(n)P(n,t)}, t>0 (2) 

being E the step operator: E k [f(n)] = f(n + k). 

Our main aim in this paper is to consider that the creation of an X particle takes a 
finite amount of time r. More specifically, we consider that the creation is a stochastic 
process initiated at a rate C(n) but, once initiated, it takes a finite amount of time r to be 
completed. Schematically: 

c D 
0^X*, X*^X, X^0. (3) 

T 

Here, X* is considered to signal the beginning of the process that, after a time r, will 
lead to X. The creation of X* is a stochastic process, but the step leading from X* to X is 
deterministic, requiring a constant time r for completion (this is indicated by a double arrow, 
while a single arrow denotes an stochastic event). We consider that the stochastic variable n 



takes into account just the number of X molecules (not including X*). Similar, althoug 



i not 



identical, processes have been considered before in the context of protein synthesis [12|, [13]. 
In a very simple manner, we can think that X* indicates the beginning of the transcription 
process of a protein from a gene, but that once initiated the transcription plus translation 
steps take a time r to be completed. In this case, the creation rate C(n) depends on n if 
there is auto inhibition or auto activation, leading to a negative or positive, respectively, 
feedback loop. 

Due to the presence of delay, the master equation of the process involves now the two- 
times probability distribution P(n, t; n', t') as: 



(E-l)[D(n)P(n,t)] + (E- 



J2C(n')P(n',t-r;n,t) 



(4) 



in'=0 

valid for t > 0. The creation term takes into consideration that the probability that a 
particle is created at time t is the sum of all contributions in which there were n' particles at 
time t — t and an X* particle was created at that time, at a rate C(n'), leading necessarily 
at time t to an X particle. This equation is the basis of our subsequent analysis. We will 
eventually consider, for the sake of concreteness, that the annihilation of particles occurs 
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through independent events at individual rate 7 and, hence, D(n) = ^n. Formally, equation 
(]!]) can be written in the form: 

= {E - l)[D(n)P(n,t)} + (E" 1 - l)[C(n,t)P(n,t)}, (5) 

with an effective time dependent rate 

C(n,t) = (C(n'),t-T\n,t), (6) 

where we have introduced the notation for the conditional average (f(ni), ti\ri2, t 2 ) = 
/( n i)-P( n i> ti\ n 2, h)- Other, non-conditional averages, will be denoted as (f(n)) t = 
f{ n )P{ n i t)- The conditional average (n,t\k,to) satisfies the evolution equation, valid 
for t > 0: 

d{n, t\k, to) 



dt 



-y(n,t\k,t ) + (C{n),t-T\k,t Q ), (7) 



with initial condition (n,to\k,to) = k. Higher order averages obey a hierarchy of equations 
which we do not need to write down for the purposes of this paper. The resolution of this 
hierarchy would allow one to compute the average value (f(n)) t of any function f(n) which 
can be expanded Taylor series of n. 

We will be mostly interested in the steady-state, where the averages (f(n)) st = 
lim t ^ l00 (/(n)) t are time independent and the conditional averages depend only on the 
time difference, (f(rii), t\n2) s t = hnit'-5>oo(/(^i), t + t'\n 2 , t'). They can be computed, re- 
spectively, from the steady-state probability distributions P s t{n) = lim^oo P(n, t) and 
P s tini 1 t\n2) = lhm/^oo P st {rii, t' + t\n 2 ,t'). Formally, the knowledge of the steady value 
C st {n) = lim^ 00 (C(n / ), t — r\n,t)) = (C(n'), — r\n) s t, allows the calculation of the steady- 



state probabilities P st (n), after imposing 8P ^"'^ = in Eq.([n]), as 



at 



Q: 



p ^) = p -(°mmfr ) = ^U^), (8) 

k=0 V ; ' k=0 

P s t{0) is fixed by the normalization condition. In the following subsections, we will discuss 
two methods to obtain the conditional averages needed for the calculation of C st {n). This 
will allow us to obtain the steady state probabilities as well as the mean value (n) st , variance 
<j 2 st and correlations Kit) = (n(n' ,t\n) st ) s t — ( n )1f 
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A. The independent-times approximation 



The first method assumes that the conditional average values do not depend on pre- 
vious history or, equivalently, that the two-times probability distribution factorizes as 
P(ni,ti,n2,t 2 ) = P(rii,ti)P(n2,t 2 ). This implies that in Eq.flB]) we can set C(n,t) = 
(C(n'))t- T , independent on n. On empirical grounds, it is expected that this approxi- 
mation will be valid for large r where the events at t and t — t can be considered inde- 
pendent, although we will show later that this is not the case. In the steady state, this 
assumption implies C(n,t) = (C(n)) st , a constant. Replacing this result in Eq.flH]) we ob- 
tain that the steady-state follows a Poisson distribution P s t(n) = e~ x ^r with x — ^ c ^ 3t . 
The, yet unknown, steady state average value is obtained through the consistency relation 
(C(n)) st = C{n)P s t{n). Once this equation is solved, the mean and variance of the 
distribution follow: (n) st = a 2 st = M^lk. 

The consistency relation can be explicitly solved in the linear case C (n) = c — en with 
the result (n) st = In the literature, and in the field of protein transcription, it is 

usually considered a negative feedback loop where the creation rate is a decreasing, non 
linear function, for example C(n) = jf^- This corresponds to a gene repressed directly by 
the protein it encodes for, in the limit where the binding and unbinding of this protein is 
fast compared to the rest of time scales of the system. The approximation is good if the 
unbinding rate of the protein from the promotor is much greater (or the order of ten times) 
than the degradation rate of the protein [l^]. In this case, the consistency relation reduces 



to: 

(n) st = - 
7 



POO 

/ dxexp[-x + (n) st (e- ex -l)l (9) 
Jo 

which, in general, needs to be solved numerically. In the limit e — > we can expand 
e~ ex — 1 = —ex to derive (n) st = t , the mean- field result. Since at the steady state 

the effective creation rate C(k,t) = (C(n)) st is constant, the process is a simple birth-death 
process and the correlations decay exponentially as K(t) = cr^e -7 '*'. 

Within this independent-times approximation, the steady state average value {n) st and 
variance a 2 st are equal (Poisson distribution) and do not depend on the delay time r. This 
is a general result that does not depend on the specific functional form for the creation rate 
C(n). As discussed before, this is naively expected to hold in the case of a large delay r. In 
the numerical simulations, however, it is observed that the fluctuations are sub-Poissonian 
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for small r and super- Poissonian for large r. The details of the simulations for this stochastic 
process including delay are given in the Appendix 1. Note that the case r = can be solved 
(exactly) by a variety of methods. Within our treatment, and according to Eq. (f5|) , for r = 
the conditional probability is C(n,t) = (C(n ; ),t\n,t) = C(n), which leads to a steady state 



distribution P st (n 



Pst(0) 

7 n ra! 



n 



n-l 
k=0 



C(k). In the non-linear case C(n) 



1+en 



this leads to: 



Pst(n) 



(n> 



Ji_ 1 (2u)n!r(n+f; 
h(2v) 



St 



a 



St 



e 

(n) st - v 2 



h(2v) 
h\{2v) 



h +1 (2v) 



(10) 
(11) 
(12) 



where v = It is possible to show that a 2 st < (n) st , a sub- Poissonian distribution in this 

case of r = 0. In the next section, we will introduce an approximation that will allow us 
to explain that fluctuations can be amplified and become super-Poissonian when we include 
time-delay terms in the process. 



B. The time-reversal invariance assumption 

One of the difficulties for the calculation of C st {n) = lim t __ !>00 (C(n / ), t — r\n,t)) = 
(C(n'),-T\n) st is that it is a correlation backwards in time, whereas Eqs. (j3J) and (JTJ) are 
only valid for t > 0. The approximation we propose in this subsection is to use a time- 
reversal invariance assumption in the steady state, namely (C(n'), — r\n) s t = {C{n'),r\n) s t. 
A simple algebra shows that a sufficient condition for this time-reversal invariance to hold 
is that the stationary probabilities satisfy: P s t(n' ,t\n)P st (n) = P st (n,t\n')P s t(n'), valid for 
all t > 0. This relation is correct in the limit t — > 0, as P st (n',dt\n) = w(n — > n')dt, the 
rate of going from n to n' particles during time dt, and it then becomes w{n — > n')P st (n) = 
w(n' — >• n)P st (n'), the detailed balance condition, which is valid for any one-step process, 
as can be derived from the master equationpj. If the process were Markovian, the detailed 
balance condition would imply the time-reversal invariance for arbitrary, finite, time t. As 
the presence of delay makes the process not Markovian, the time-reversal invariance is an 
assumption whose validity and implications need to be checked. In Fig. ([1]) we plot the 
correlations (n,r\k) st and (n,—r\k) st as a function of k, using a negative feedback loop 
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C{ n ) = f° r two different sets of parameters. In the same figure we plot the stationary 
probability distribution P st (k). As it can be seen from this figure, it is not true that these 
two correlations are identical for all values of k. However, it has to be noticed that the larger 
discrepancies occur for those values of k which have a low probability of appearance. 




FIG. 1: Conditional averages in the steady state, (n, r\k) s t (+ symbols) and (n, — T\k) s t (x symbols) 

coming from numerical simulations of the process with delay schematized in Eq. (|3j) using a creation 

rate C(n) = , with t = 10, eo = 1, cq = 3 and two different values of Q = 50 (top) and = 5 

i+ n n 

(bottom). In the same figures, we also plot with □ symbols, the (arbitrarily rescaled) stationary 
probability distribution P st (k). Note that the discrepancy between {n,r\k) st and {n, — T\k) st is 
larger in those cases that the particular value of k is less probable. 

Once this time reversal invariance assumption has been adopted, to compute C s t(n) = 
{C(n'),T\n) st , one could solve the hierarchy of equations for the moments with the appro- 
priate initial condition. This could be done, for instance, if the creation rate C(n) were a 
linear function C(n) = c — en. However, as discussed before, most interesting cases consider 
a negative feedback with a nondinear rate C(n). In this case, one can not, in general, close 
that hierarchy of equations and one needs approximate methods to find C st {n), such as, 
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for example, the Gaussian closure 15[ . In the following, and in the spirit of van Kampen's 
expansion 1,| , we will linearize the equations assuming that the variable n has a deterministic 
contribution of order Q (a large parameter of the system, typically the system volume) and 
a fluctuating part of order i.e. n = f20 + fl^^. Although it is possible to deal with 
the most general case, we will restrict ourselves to the case where the creation rate satisfies 
the following scaling with system size C{n) = VL<^> (^), so one can expand C{n) around the 
macroscopic component: 

C(n) = n ($ (0) + fi- 1/2 ^'(0)e + •••), (13) 

so that 

(C(n'),t'\n, t) = n* (0(0) + ^M^OXC, t% t), (14) 

with f = Q-^n - Q 1 / 2 ^. 

We replace ansatz ( |T3l) in the evolution for the first moment (171), and equate the powers 
of Q to find that the deterministic (macroscopic) and stochastic contributions to n satisfy: 

d<f>(t) 



dt 

d(e,t%t) 



- 7 0(t)+^(0(t-r)), (15) 
-^',t%t)+^(<f>(t-T)){e,t'-T\U) (16) 



dt' 

Equation (j!5p for the macroscopic component is, in general, a nonlinear delayed differential 
equation which might be difficult to solve. However, the steady state value 4> st is readily 
accessible as the solution of 70 st = ${<p st ). The stability of this fixed point is found by 
linearization around it. A standard analysis of the resulting linear delay differential equation, 
tells us that a sufficient (but not necessary) condition for stability is \a\ < 7, where we have 
defined a = —<&' (0 s t). 

Once in the steady state, we replace (pit) by its stationary value 4> st and Eq. ffTBT) becomes 
a delay linear differential equation with constant coefficients, and we are looking for the 
time-symmetric solution of this equation satisfying the initial condition (£', t\£,t) = £. This 
can be written as (£',£ + A|£, t) = £/(A), being f(t) the symmetric solution f(—t) = 
fit) of the equation /(£) = —jf(t) — af(t — r) and /(0) = 1 (see Appendix 2). From 
Eq.(Hl]) we get the effective creation rate C{n) = VL<P{<p st ) + fi 1/2 ^ / (0^)^/( r ) = ^4>st(l ~ 
®'(<t>st)f(r)) + <P'((j) s t)f(T)n after replacing f = tl^n - fi 1 / 2 ^ and 70 st = <P((j) st ). From 
Eq.® one can obtain the steady-state probabilities P s t(n). Their functional form depends 



8 



on the sign of &(4> at )f(T): (i) If ^'(0 si )/(r) < 0, the distribution is a binomial distribution 

p st(n) = Ov n {i - p) M ~ n with p = ;%t)\ T (l) and M = m ^ { -nl)f(r) ~ x ) and 

< n < M; (ii) if <l?\(j) st ) f (t) = 0, the distribution has a Poisson form P s t(n) = e _x ^j- 
with x — (hi) finally, if <P'((f) s t)f(r) > 0, the distribution is a negative binomial, 

Pst{n) = r+r'X 1 - <?) A V, with q = and M = Q<P(4> st ) [ W ^ m - l). In 

all cases, however, they can be approximated up to terms of order fi -1 / 2 by a Gaussian 
distribution. Despite the differences in the functional form, in all three cases the mean value 
and variance are given by: 



[n) st = (17) 

(18) 



2 ( n )st 



st l- 7 -^'(0 st )/(r)' 
An equivalent expression for the variance taking as a starting point a linear Langevin dif- 
ferential equation including delay was obtained in jsl, Q. 

In the case of a negative feedback loop, it is a = —$'((j) s t) > 0. It can then be seen from 
the expression in the Appendix 2 that /(r) monotonically decreases from the value 1 at 
t = to the value — < at r — > oo (A = a/7 2 — a 2 , see Appendix 2, recall that |a| < 7 
is a sufficient condition for the stability of the fixed point <p s t)- In this case the fluctuations 
are sub-Poissonian if f(r) > (small r) and super-Poissonian if /(r) < (large r). The 
threshold between the two cases is the value Tp at which f(rp) = or Tp = —A -1 ln£ in the 
notation of the Appendix 2. As explained before, the probability distribution is binomial 
for r < Tp, Poissonian for r = Tp and a negative binomial for r > Tp. 

In the case of positive feedback, a = —<P'((j) s t) < 0, /(r) monotonically decreases from 
1 at r = to > at r — y 00, and in this case the fluctuations are always super- 

Poissonian, but their magnitude is reduced as the delay is increased. The steady-state 
probability distribution is always a negative binomial distribution. 

We conclude that the delay can have opposite effects: in a negative feedback loop it 
enhances the fluctuations, whereas in a positive feedback loop it reduces them. On the 
other hand, it is well known that, in the non-delay scenario, a negative feedback reduces 

n 

the magnitude of the fluctuations [17J when compared to the n-independent creation rate. 
We find it remarkable that the presence of delay can reverse the usual fluctuations-reducing 
effect of the negative feedback loop, and, instead, enhance the fluctuations. 

The correlations in the steady state can be obtained from K{t) = (n{n',t\n) st ) st — (n) 2 sV 



as: 

2 



K(t) = a 2 st f(t). (19) 

Note that, as can be seen from the alternative definition K(t) = ]im t i^^(n(t+t')n(t'))—(n)^ t , 
the correlation function is a time-symmetric function K(—t) = K(t). However, and contrary 
to previous assumptions [3], this does not imply that the conditional expectation value 
(n',t\n)st has to be a symmetric function. In fact, it is not for an arbitrary value of n, as 
shown in FigJT] 

We apply these results to specific functional dependences of C(n). Let us first comment 
that in the linear case C(n) = c — en, Eq.([7|) is already a closed equation and our treatment, 
not surprisingly, can be carried out without assuming the expansion f[T^]) . However, we do 
not find this case very interesting as it turns out that the problem is ill-defined as the rate 
C(n) might become negative when the number of molecules n exceeds c/e. 



18| , is the rate C ( 



A more interesting case, used in the protein transcription problem|lN|. is t lie rate ( '( // ) 
that we write in the form C{n) = VL<£> (^) with $(z) = and c = c/Q, e = eQ 

where Q is a large parameter, typically proportional to the cell volume. This corresponds 
to a negative feedback loop. Note that the condition |a| < 7 is always satisfied for such a 
creation rate and the steady state 4> st is always stable no matter how large the delay time r. 

In Fig.(T5|) we compare the average and variance obtained from numerical simulations 
with those obtained from the theoretical analysis. The agreement is, in general, very good 
and improves as Q becomes large. In Fig. ([3]) we compare the correlation function obtained 
numerically with the analytical expression [19j Its non-monotonic character due to the delay 
is apparent. The value of the correlation at t = r is not negligible, compromising the validity 
of the independent-times approximation. 

For C(n) = — c °" n with I > 1, (negative feedback loop with cooperativity) the equation 

1+e o(f!) 



for the macroscopic variable ([15]) has a Hopf bifurcation into a limit cycle attractor. For 
parameters below the Hopf bifurcation, the situation is qualitatively equal to the previous 
case, and the discussion applies. For parameters above the Hopf bifurcation, the system 
becomes oscillatory so the assumption of steady state is not valid, and the results obtained 
here are not directly applicable. 
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FIG. 2: Steady state average (n) st (dashed lines) and variance a 2 st (full lines), for process denned 
in Q, as a function of the delay time r, for a creation rate C(n) = f°ijj with Co = 3 (upper 
part of each panel) and cq = 1 (lower part of each panel), and two system sizes (Q) (upper and 
lower panel) and cq = 1 in both cases. In each case, we plot with symbols the results coming from 
numerical simulations and by lines the theoretical expressions, Eqs. (|17p and ()18p . 

III. DISTRIBUTED DELAY 

In general terms, it is more realistic to consider that the delay that each individual event 
takes to be completed is a fluctuating quantity following some probability distribution rather 
than taking a fixed value. This is definitely the case in genetic networks, where transcription 
and translation times can be broadly distributed [19]. In this section we exemplify how to 
apply the method developed before in the case of a stochastic delayed production process 
including distributed delay. 

We consider again the process schematized in ([3]), but now we consider that the delay 
time t is a stochastic variable with some probability distribution p(r). For simplicity, we 
consider that the delay times for all individual reactions are independent and identically 
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FIG. 3: Correlation function in the steady state, for the delayed process ([3|) with creation rate 



C{n) 



Simulations (circles) and theory, equation (fT9|) (solid line). 



distributed. Now, the master equation for the process is: 
dP(n,t) 



dt 



(E-l)[D{n)P(n,t)} + (E- l -l) 



/ drp(r)C(n')P(n',t-r;n,t) 
,n'=0 J 



(E - l)[D(n)P(n,t)] + (E~ - l)[C(n,t)P(n,t)] 



(20) 



with C(n,t) = J dTp(r)(C(n'),t — r\n,t) = f drp(r)C(ri,t;T). Hence, it is possible to follow 
formally the method of last section simply replacing C by C and we skip the details of the 
calculation. The mean value, variance and correlation function are given by: 



(n) st 

2 



00. 



st ) 



K(t) 



l-J-^'(4> st )JdTp(T)f(Ty 



being f(t) the solution of the integro-differential equation 

df{t) 



dt 



-!f{t)+& (<P st ) J drp(r)f(t-T) 



(21) 
(22) 

(23) 
(24) 



satisfying f(-t) = f(t) and /(0) = 1. 

There is no general method that can be applied to find the solution of this compli- 
cated equation. A reduction to a set of linear differential equations can be achieved if 
we adopt the Gamma probability distribution: p(r; k) = Ar k ~ x e~^' T ', depending on two 
parameters: k and r. The average value is r and the root-mean-square is a T = In- 
creasing k for fixed f decreases the fluctuations of r, and in the limit k — > oo the distri- 
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bution approaches a Dirac-delta and r becomes a deterministic variable (fixed delay, cor- 
responding to the case analyzed in the previous section). The alternative solution method, 



known as the linear-chain trick [20|, begins by defining a family of time-dependent functions 
Zi(t) = J drp(r; l)f(t — T),l — l,...,k. After some algebra, one can prove that f[2~21) is 
equivalent to the system of linear ordinary differential equations: 

df{t) 



- 7 /(t)+<2>'(^ t )Z fe (t), (25) 
(/(*)- Zi), (26) 



dt 

dZ\ k 
~d~t ~ T' 

^ = -(Zi-i-Z^, l = 2,...,k. (27) 
at r 

which, besides f(0) = 1, require a set of initial conditions for Z\{t = 0), I = 1, . . . , k. These 
can be determined in a self-consistent manner. First, note that the symmetry condition 
f(t) = f(-t) implies: 

Zi{t = 0) = j drp(r; l)f(r), l = l,...,k. (28) 

One then solves ( 1251127]) with arbitrary initial conditions for Z\(t = 0) and imposes fl28|) . 
This yields an algebraic system of k linear equations for Z\{t = 0). The solution of the 
linear differential equations ( 12511271) and the solution of the algebraic equations ( 1281) can be 
obtained, either analytically for small k, or numerically, but with a very high precision, 
for large k. Note that in order to compute the variance, Eq. (1221) . all we need to know is 
JdT P (r;k)f(r) = Z k (t = 0). 

In Fig.(jl]) we plot the ratio a 2 st / (n) st as a function of a T for fixed mean delay r. We see 
that as the delay distribution becomes wider (decreasing k), the fluctuations of the process 
decrease, so that the effect of the delay becomes less important. The results for the Gamma 
probability distribution are qualitatively equal to other distributions for the delay times such 
as uniform or Gaussian (truncated in order not to produce negative values). This ressults 
suggests that a natural or artificial system should have a rather precise delay if it is to make 
use of the effects that delay induces in the fluctuations, or it should have an irregular delay 
to avoid those effects. 
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FIG. 4: Variance normalized to the mean value o%/(n)st) for the process with distributed delay 

defined in ([3]) and a creation rate C(n) = 1 ^°j? n , for a delay distributed according to a gamma 

distribution p(r; k) = Ar k ~ 1 e~^ T , as a function of the relative size of the fluctuations in the delay 

-L = Ar 1 / 2 . Results coming from numerical simulations (o) and from the theoretical method (x) 
r 

as explained in the main text. 



IV. TRANSCRIPTION-TRANSLATION MODEL 



So far, we have consider simple one-step birth and death processes. In the context of 
gene regulation, however, the protein production involves two major steps (transcription 
and translation) and it is well known that the combined effect of the two steps can enhance 

n 

significantly protein fluctuations [17j. In this section we study the effect of delay in a more 
elaborated model for protein levels than the one considered previously, including explicitly 
the transcription (creation of mRNA from the DNA) and translation (creation of the protein 
from the mRNA) steps. The process can be schematized as follows: 

C 0J 7n 7, 7 i 

0^r*^r, y^x*^x, x^0, v^0. (29) 

Now X corresponds to the protein (with n the current number) and Y to the mRNA. We 
denote by m the number of mRNA molecules at time t — t%. In doing so, the translational 
delays T\ and t 2 can be absorbed in a total delay r = T\ + t 2 . The master equation for the 
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process is: 

dP(m, n, t) 
dt 



(E n - 1) [y n nP(m, n, t)] + (E m - 1) [ lm mP(m, n, t)] (30) 

oo 

C(n')P(n', t — r;m, n, t) 



+ (E- 1 - 1) [ujmP(n, m, t)] + {E' 1 - I) 

_n'=0 

being E n and E m the step operators for the number of proteins, n, and the number of 
mRNA, m, respectively. As before, we will allow for feedback loops by letting the creation 
rate C to become a function on n. For simplicity, though, the translation rate u, as well as 
the degradations rates 7„ and 7 m will be considered constant. 

The general formal expression for the stationary solution of the master equation ( 1301) is 
not known. To proceed in this case, we will apply van Kampen's expansion, which assumes 
both n and m to be split in deterministic and stochastic contributions as n = 
and m = Vt<p m + Vl}/ 2 ^ m . The probability density function n(£ n ,£ m ) for the stochastic 
variables satisfies a Fokker-Planck equation that is found by expanding the master equation 
in powers of Q: 

dmm d f n ^ = ^{h m u-f'(ut-rme n ,t-T\u^n,t)}u} 

1 c) 

+ ■= hm<f>m + f{<t>n{t - t))] ^11 + -^{[jntn ~ wfm] n} 
1 d 2 

+ ~ [ln<Pn + U<Pm] ^2 II. (31) 

The deterministic contributions <p n , (fi m and the averages of the fluctuation terms obey the 
following system of delayed differential equations: 

-7m0m + $(0n(*-T)), (32) 



dt 



dt 

dt' 



r = ~7n0n + W0 m , (33) 
-lm(C^n,Ut) + $'(0n(t " r)){^t' - T^ n ^ m ,t), (34) 
-ln(&,t\€n,€m,t) + Uj{£' m ,t'\£ n} £ m ,t) . (35) 



dt' 

The solutions for the average of the fluctuations with appropriate initial conditions, after 
replacing (j) m (t) and <f> n (t) by their stationary values 0„ iSi and 4> m ,st coming from the fixed- 
point solution of Eqs. ( 13211331 can be solved under the assumption of time- reversal invariance, 
to obtain: 

(C t\U, Dst = /„(*)£„ + f m {t)im (36) 
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(see Appendix 2 for explicit expressions of the functions f n (t) and f m (t)). We replace again 
4> m (t) and 4> n (t) by 4>n,st an d 4> m ,st and use the time reversal approximation (£' n , — r|£ m , £ n )st = 
(Cn> r l£m> ^n)st to reduce Eq. (l3"Tj) to a linear Fokker-Planck equation whose solution is well 
known to be a Gaussian distribution[l]. The corresponding steady state values for the 
average and fluctuations in protein levels are given by: 



n )st — ^</>rc,st 
2 

1 + 



n,st 



7m 



l-yn(r)(l + ^/ m (r)) 



<n> 



1 + i + ^(-)l + ^(^(r) + / m 



(37) 
(38) 



In the case of no delay (r = 0), this expression reduces to the one obtained in [17J. In Fig. dSJ) 
we compare the average and variance of this transcription-translation model as a function 
of the delay for a creation rate of the form C(n) = 1 '^^_ n ■ Again, in this negative feedback 
loop setting, the delay significantly enhances the fluctuations, up to a level well over the 
value without feedback (marked in the figure by a dashed line), leaving the mean value (n) st 
essentially unchanged. So again in this case, the delay reverts the effect of the negative 
feedback, from fluctuation-reducing (for low values of the delay) to fluctuation-amplifying 
(for large values of the delay). 



£1=50, £ =1, co=30, y =1, y =5, c =3 

(1 'n 'm 




FIG. 5: Stationary values for the average (n) s t and variance a st for the protein levels as a function 
of the total delay, for the transcription-translation model schematized in (|29p for a creation rate 



of the form C{ri) 



Values from numerical simulations (symbols) and theory (solid lines, 



Eq. (I38p ). Values of parameters in top of figure. The dashed line corresponds to the variance of a 
system without feedback, with the same average. 
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V. DISCUSSION 



We have studied stochastic processes with discrete variables in continuous time that 
include delay. We have shown that the combined effect of feedback and delay gives rise 
to nontrivial results. When a stochastic process has negative feedback, the fluctuations 
are decreased; if this feedback is delayed, the fluctuations can be are actually enhanced, 
depending on the magnitude of the delay. A positive feedback loop enhances the fluctuations, 
but if the feedback is delayed, this enhancement is decreased. We have also shown that this 
effect of the delay is less apparent if the delay itself has relative large fluctuations, so for 
this mechanism to work, the delay has to be controlled precisely. This may be relevant for 
example in gene-regulatory networks, where delay times are typically broadly distributed 
but several regulatory mechanisms may act to control this 19j. The analytical theory allows 



us to understand and predict this phenomenology in a general way. We have also shown 
that the assumption of de-correlation of times t and t — r for large delays is not justified a 
priori, since the correlation function is typically non-monotonically decreasing, with peaks 
at multiples of the delay. Finally, we have pointed out that systems with delay are not, in 
general, statistically invariant under time reversal over the steady state, even if they fulfill 
the detailed balance condition. 



Appendix 1: numerical simulations 

realizations of the process, we use the following modification of the 



To perform numerica 
Gillespie algorithm 21 



22|: 



1: Initialize the state of the system, setting, e.g. n — 0. 

2: Compute the reaction rates C(n) and jn. Obtain a number At exponentially dis- 
tributed with average l/(C(n) + jn). 

3: If t + At is larger than the time of the next scheduled delayed reaction, go to step 
4. Otherwise, update time from t to t + At and obtain which kind of process (creation or 
degradation) will take place. To do so, generate a uniform random number between and 
1. If this number is smaller than jn/(C(n) + jn), set n — > n — 1; otherwise add an entry in 
the list of scheduled creation processes to happen at time t + r. Go to step 2. 

4: Update the time to that of the next scheduled reaction. Set n — y n+ 1. Go to step 2. 
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This procedure is statistically exact, as the original Gillespie algorithm in the case of 
non-delayed reactions. 

In the case with delay, the time until the next reaction is exponentially distributed, with 
average C(n) + 771, only if the state of the system doesn't change during this interval (due 
to a scheduled delayed reaction). This happens with probability 1 — e~^ c ^ +in ^ tr (with 
t + t T the time of the next scheduled delayed reaction). The algorithm fulfills this, since the 
probability that step 3 is completed is precisely 1 — e -{ c i n )+^ n ) t r _ Once a reaction has taken 
place (delayed or not) the time for the next reaction is again exponentially distributed as 
long as no delayed reaction takes place, and the procedure can be iterated. 



Appendix 2: solution of the delay-linear equations 



We consider the following linear delayed differential equation: 

df(t) 



dt 



-af(t-r)-jf(t). 



(39) 



We are looking for a symmetric solution /(— £) = f(t). We summarize here for completeness 



the treatment of reference 12]. We make the ansatz f(t) = ae A '*' + be A '*', valid only for 
—t < t < r. Inserting in ( 1391) . equating the coefficients of e xt and e~ xt , and imposing 



/(0) = 1, we obtain X,a,b. Once we know f(t) for \t\ < r, we can obtain f(t) for \t\ > r 
iteratively integrating (139]) . The solution for t > is: 

7 — A 



A = a/7 2 — a 2 



c 



a 



f(t) 



e —\t_£ e \(t—T) 

l-Ce" AT ' 



if < t < t 



(40) 



-7(t— fcr) 



/(hj-a/^i'f-rje^), if k r < t < (k + l)r, k = 1, 2, 



Note that /(t) 



1-Ce- AT ' 



Using the symbolic manipulation program Mathematica 



23| to 



perform the integrals of the iterative process, we have been able to find explicit expressions 
for f(t) up to \t\ < lOr. 

We apply a similar approach to the case of two coupled linear delayed differential equa- 



ls 



tions: 



dx m (t) 

dt 
dx n jt) 

dt 



1m%m{t) O^X m (t 7"), 

(t) +wx m {t). 



Due to the linearity, the solution has the form: 



X n (t) = X n (Q)f n (t) + X m (0)f m (t), 



(41) 
(42) 

(43) 



with / n (0) = 1, / m (0) = 0. To find this solution, we use the ansatz x n (t) = aie A+ '*' + 
5 ie -A+|*| + a2 e A -l*l + 6 2 e A - |t| , x m (t) = cie A +l*l + c 2 e~ x + w + rf 1 e A - |t| + d 2 e x -^, for -r < i < r. 
Equating the coefficients of the exponentials and imposing the initial condition we obtain 
the expression valid in < t < r: 

1 



fn(t) = 



In 



b-(t) +x l + b-(t) 



In- 



bit) 

i -Kit) 



fm(t) 



CO 



b(t) 
l-b + (t) 



b(t) 

bit) 



(e x+t -b + it)e~ x+t ) 



(e A -*-6_(t)e' A -*) , 



bit) 



[e x -t - b_it)e- x -t] - b X+t ~ b + (t)e~^] , 



A 



± 



lm _ ln ± -,/7 
2 2 v ' 



6 ± (0 
6(0 



A| + (7 m + 7n)A ± + 7^ 



wa 



A_(l + 6_(t))(l - &+(*)) - A+(l - b-(t))(l + b+(t)). 



(44) 
(45) 
(46) 

(47) 

(48) 
(49) 
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